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The iterative phase retrieval problem for complex-valued objects from Fourier transform magni¬ 
tude data is known to suffer from the twin image problem. In particular, when the object support 
is centro-symmetric, the iterative solution often stagnates such that the resultant complex image 
contains the features of both the desired solution and its inverted and complex-conjugated replica. 
The conventional approach to address the twin image problem is to modify the object support dur¬ 
ing initial iterations which can possibly lead to elimination of one of the twin images. However, at 
present there seems to be no deterministic procedure to make sure that the twin image will always 
be very weak or absent. In this work we make an important observation that the ideal solution 
without the twin image is typically more sparse (in some suitable transform domain) as compared 
to the stagnated solution containing the twin image. We further show that introducing a sparsity 
enhancing step in the iterative algorithm can address the twin image problem without the need to 
change the object support throughout the iterative process even when the object support is centro- 
symmetric. In a simulation study, we use binary and gray-scale pure phase objects and illustrate the 
effectiveness of the sparsity assisted phase recovery in the context of the twin image problem. The 
results have important implications for a wide range of topics in Physics where the phase retrieval 
problem plays a central role. 
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The problem of phase measurement carries enormous 
importance in many disciplines of Physics such as astro¬ 
nomical imaging [1], microscopy of transparent biologi¬ 
cal cells [2], coherent X-ray imaging lam, electron mi¬ 
croscopy [ 7 ], and imaging using ultrashort X-ray pulses 
[8], to name a few applications. In the context of coherent 
scattering of light waves, it is known for long time that 
the phase variations carry much more useful information 
as compared to the amplitude variations [9]. Phase is 
not directly measurable but can only be inferred compu¬ 
tationally using interferometric methods or by applying 
iterative phase retrieval algorithms to the direct far-field 
diffraction intensity measurements. The iterative phase 
retrieval based techniques using some prior information 
about the object to be imaged are attractive as they do 
not require complicated interferometer setups that may 
be difficult to build when using X-rays or electron beams. 
Following the early work by Gerchberg and Saxton m, 
and Fienup mm, iterative phase retrieval remains an 
active research area m- One of the most difficult prob¬ 
lems associated with iterative phase retrieval is that of 
recovering a complex valued object g{x, y) from the mag¬ 
nitude \G{fx,fy)\ of its 2D Fourier transform ITS]. 
The amplitude pattern \G{fx^fy)\ in this case does not 
possess any symmetry. Further, since the two functions 
g{x^ y) and g'^{—x, —y) have the same Fourier magnitude, 
iterative phase retrieval methods are known to suffer from 
the twin image problem m, where the recovered solu¬ 
tion consists of features of both g{x, y) and the twin im¬ 
age ^*(—X, —^). Although a linear combination of the 
form tg{x, y) {1 — t)^*(—x, —y) with t G (0,1) does not 
have Fourier transform magnitude equal to \G{fx,fy)\^ 
the iterative solution is known to stagnate and not make 
progress towards the desired solution containing only one 


of the twin images corresponding to t = 0 or t = I. The 
problem gets much more severe when the object support 
is centro-symmetric which is often a natural choice for 
the support constraint to be used along with experimen¬ 
tal Fourier magnitude data. Our aim in this Letter is to 
provide an image sparsity based method for addressing 
the twin image problem. We note that image sparsity 
provides an independent criterion other than the usual 
support constraint which can help eliminate the twin im¬ 
age from the iterative reconstruction. 

A few solutions have been proposed in the literature 
to address the twin image problem. The most promi¬ 
nent approach to avoid the twin image is to truncate 
the object support to a non-centrosymmetric window for 
a few initial phase retrieval iterations [16]. With this 
procedure, there is a good chance that one of the two 
possible solutions gains prominence and using the full 
support in the subsequent iterations does not bring the 
twin image back. It has also been shown m that if the 
Fourier transform is sufficiently over-sampled, then using 
the difference map algorithm [18] with an appropriate pa¬ 
rameter choice, the twin image problem can be avoided. 
A moving overlapping aperture approach where multi¬ 
ple intensity measurements are made [19] has also been 
suggested for eliminating the twin image. More recently 
an interesting observation has been made [20] that at the 
stagnation stage, the retrieved Fourier domain phase gets 
divided into two regions. One of the regions corresponds 
to the upright image and the other region corresponds 
to its twin thus to some extent justifying the truncated 
support based solution. We believe that while the above 
solutions do offer some insights into the twin image prob¬ 
lem, there seems to be no deterministic procedure that 
can drive the solution out of the stagnation stage when 
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the Fourier magnitude data has been sampled just ade¬ 
quately as per the Nyquist criterion. In this Letter we 
present a solution to the twin image problem without the 
need to change the support constraint during iterations 
or the use of oversampling in the Fourier transform space. 
As shown here when an image sparsity enhancing step is 
added to a standard phase retrieval approach such as the 
hybrid input-output method (HIO) [12], the twin image 
always appears to get eliminated (or possibly weakened) 
even when the object support is centro-symmetric. While 
sparsity based algorithms for phase retrieval have been 
proposed in recent years, they do not specifically address 
the twin image problem so far [21]. 

We start by making an important observation that the 
ideal solution g{x^y) or g*{—x, —y) is typically sparse in 
some suitable transform domain as compared to the stag¬ 
nated solution containing the twin image. As an illustra¬ 
tion we use a unit amplitude pure phase object g{x, y) 
whose amplitude and phase are shown in Fig. 1 (a), (b) 
respectively. A phase step of 27r/3 is used in the binary 
phase object. The corresponding Fourier transform mag¬ 
nitude data \G{fx^ fy)\ is shown in Fig. 1 (c). The image 
size used is 512 x 512 and the object support is assumed 
to be the central 250 x 250 pixels which is slightly smaller 
than half the image size so as to ensure adequate sam¬ 
pling in Fourier transform space as per the Nyquist crite¬ 
rion. Further, the support window has been purposefully 
made centro-symmetric. A typical stagnated phase so¬ 
lution containing the twin image is shown in Fig. 1(d). 
Here we used 500 iterations of the HIO method for phase 
retrieval with object support as the only constraint. In 
this and further illustrations in this paper, we observed 
that by the end of 500 iterations the phase solution was 
stabilized. For completeness we mention that in our nu¬ 
merical study for the binary phase object above, truncat¬ 
ing the support to a right triangular window bounded by 
two sides and diagonal of the square support window for 
a fixed number (= 10) of initial iterations did not elim¬ 
inate the twin image in a series of 20 runs of the HIO 
method with different initial guess for the phase of the 
Fourier transform G{fx,fy)- The initial guesses for the 
phase were uniform random phase patterns with phase 
values distributed in [0,27r].The truncation of the sup¬ 
port window during initial iterations thus does not seem 
to guarantee the twin image removal and a more deter¬ 
ministic approach is required to address this problem. 

A suitable sparsity measure for the binary phase object 
in this example is the total variation (TV) [22] defined 
as the L 1-norm of the image gradient: 

TV{g) = \\Vg\U= Y. U) 

2=all pixels 


The numerical value of the TV of the image g{x,y) in 
Fig. l(a),(b) is seen to be 1.92 x 10^ and the TV of 
the stagnated solution is numerically equal to 8.99 x 10^. 



FIG. 1. (a) Amplitude and (b) phase of unit amplitude pure 
phase object g{x,y). Image size is 512 x 512 and the object 
support as seen in (a) is central 250 x 250 pixels. The phase 
step in the binary phase pattern in (b) is equal to 27r/3. (c) 
Fourier transform magnitude displayed as \G{fx, to 

suit the display, (d) a typical stagnated phase solution using 
500 iterations of the HIO method showing the twin image 
problem. 


In 20 consecutive runs of the HIO algorithm (with 500 
iterations each) starting with a new realization of the 
initial random phase map as a guess, we observed that 
the twin image problem always occurred for the HTD 
phase object and that the TV values of the stagnated 
solutions were clustered in a narrow range with mean 
equal to 9.01 (±0.14) x 10^. In the present case, a sparsity 
measure such as TV therefore appears to be a suitable 
criterion to distinguish a stagnated solution from a non- 
stagnated solution. Further, the above observation also 
suggests that a TV reducing step incorporated in the HIO 
algorithm can provide us a way out of the twin image 
stagnation problem as we will demonstrate next. 

The iterative algorithm used is shown schematically 
in Fig. 2 and described in detail as follows. Given the 
Fourier transform magnitude |G(/cc, /y)|, we start with a 
random phase function 0o(/cc, /?/)with pixel phase values 
uniformly distributed in [0, 27r] as a guess for the phase 
of G{fx, fy)‘ The following steps are then carried out in 
the (n ± l)-th iteration. 

1. Compute the inverse Fourier transform gn = 

=-F-M|G|exp(i0„)]. 

2. Update the solution g^ using the standard HIO 









3 


method as: 


9n+i = Qn for(a;, y) & C 

= 9n- P9n for(a:;, y) ^ C, (2) 

where C denotes the support constraint and the pa¬ 
rameter (3 is typically selected to be in the interval 
(0.5,1) (we used p = 0.9). 

3. Update Qn+i by using a fixed number Ntv of gra¬ 
dient descent steps of the following form for total 
variation reduction: 

(3) 

with Qn+i = ^n+ 1 - Here the functional gradient of 
the TV is defined as: 

V/Tn/) = -V-(^). (4) 

The step size t is determined in each iteration by 
a backtracking line search [23]. In our illustration 
we have used Ntv = 30 gradient descent steps for 
TV reduction. Also the TV reduction procedure 
was applied only to the part of the image limited 
by the support constraint rather than to the full 
image at this intermediate stage. More advanced 
recent methods for TV reduction may also be used 
instead of this step [SUES]. 

4. Set5n+i=5^-fr. 

5. Calculate the forward Fourier transform Gn+i = 
^[dn+i] and replace its magnitude with the known 
Fourier transform magnitude \G\ leaving its phase 
unchanged, so that, 

Gn+i = \G\exp[i(j)„+i], (5) 

where (/)„+i = arg(G„+i). 

In a sequence of 20 runs of the modified HIO algorithm 
above with a TV-reducing step and with the same ini¬ 
tial random guesses (l)o{fx^fy) as before, we found that 
the twin image was never perceptible and g{x,y) or 
^*(—X, —y) were reconstructed with almost equal proba¬ 
bility. Two of the typical phase reconstructions using the 
modified HIO method incorporating a TV reducing step 
are shown in Fig. 3 (a), (b) respectively where the part 
of the reconstructed image outside the support area has 
been zeroed out. The reconstructed phase maps have a 
constant phase shift relative to the ground truth phase 
map in Fig. 1(b) which is inconsequential. The numeri¬ 
cal TV values of the reconstructed complex images in the 
20 runs have a mean value of 2.02(±0.17) x 10^ which is 
close to the TV value for the ground truth image as stated 
above. While TV is a suitable sparsity criterion for the 
special case of the binary phase objects as used in the 



FIG. 2. Flowchart of the modified HIO algorithm for phase 
retrieval including TV reducing step and object support con¬ 
straint. 



FIG. 3. (a),(b) Two of the typical phase reconstructions using 
modified HIO method with TV reducing step. 


illustration in Fig.s 1 and 3, the phase recovery results 
provide us an indication that sparsity of images can in 
general be exploited for avoiding the twin image problem. 
It is now a well accepted fact that most natural images 
are sparse in some suitable transform domain. In recent 
times image sparsity has been effectively utilized for sig¬ 
nal reconstruction from “incomplete” data as in the class 
of algorithms known as compressive sensing [26]. While 
the present problem is not exactly a classic compressive 
sensing problem, we observe that image sparsity is indeed 
a useful concept here as well. 

In the remaining part of this paper we will use a gray¬ 
scale phase object so that our sparsity based ideas can 
be generalized to a wider class of applications of iterative 
phase retrieval. We use a modified Huber penalty [27]- 
[29] as the sparsity measure in this case instead of the 
TV penalty. The modified Huber penalty we use for the 
gray-scale phase object is defined as: 


H{9) = 

2=all pixels 



-1 


(6) 
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Here 6 is a tuning parameter and \\/gi\‘^ = \Vxgi\‘^ + 
For the image pixels where the gradient mag¬ 
nitude >> S‘^ the penalty is similar to the TV 

penalty function whereas when << the penalty 

function takes the form \\/gi\‘^/2S^. The quadratic 
penalty in image gradient is known to favor image 
smoothness [29]. Huber penalty is therefore a gener¬ 
alization of the TV penalty that is suitable for objects 
containing smooth regions as well as edge-like features. 
The parameter S may be selected based on the statistics 
of the gradient magnitudes over all pixels in the image. 
We selected 6 to be equal to the median of the gradient 
magnitudes over all the pixels within the support win¬ 
dow in each gradient descent step. Gradient magnitudes 
much larger than S then correspond to edges which are 
preserved and those much smaller than S correspond to 
small local oscillations that get smoothed out in the Hu¬ 
ber reducing step. 

Figure 4(a) shows phase map of the unit amplitude 
Lena phase object that we used as a gray-scale phase ob¬ 
ject. The image size is again 512 x 512 and the object 
support is taken to be the central 250 x 250 window which 
is centro-symmetric. The phase values are scaled in the 
range [0,57r/6] for this illustration. A typical stagnated 
phase recovery using 500 iterations of the HIO method 
consisting of the twin image is shown in Fig. 4(b). Fig¬ 
ures 4 (c), (d) show phase recoveries using 500 iterations 
of the modified HIO method with TV and Huber reduc¬ 
ing steps respectively. The same random phase map has 
been used as initial guess for all the three phase recoveries 
in Fig.s 4(b)-(d). As expected the TV penalty leads to 
flattening of the solution and the gray-scale features are 
better preserved when the Huber penalty is used. Once 
again in a sequence of 20 runs of the HIO method with 
different realization of 0 o(/cc 5 fy)^ the twin image was ab¬ 
sent (or very weak) for 4 of the runs and was present to a 
varying degree in rest of the 16 runs, which is consistent 
with observations in [20|. When the same (l)o{fx^fy) re¬ 
alizations were used, the twin image was however never 
perceptible with the modified HIO method including the 
TV or Huber penalty. 

In conclusion, we have demonstrated that the twin im¬ 
age problem in phase retrieval can be addressed in a de¬ 
terministic manner if the well known algorithms such as 
the HIO method are modified to include a sparsity en¬ 
hancing step. In our numerical experiments with binary 
and gray-scale phase objects, we observe that the twin 
image can be eliminated even with a centro-symmetric 
object support and when the Fourier space sampling is 
just adequate as per the Nyquist criterion. While we have 
used the TV and the modified Huber penalty as sparsity 
measures for illustration, we believe that there is a vari¬ 
ety of other sparsity measures such as the wavelet domain 
sparsity that may also be utilized effectively for address¬ 
ing the twin image problem. The specific choice of the 
sparsity measure will always depend on the problem or 



FIG. 4. (a) Phase map of the unit amplitude Lena phase 

object, phase values scaled to the range [0,57r/6], (b) phase 
recovery using 500 iterations of HIO method, (c) phase recov¬ 
ery with modified HIO method with TV penalty, (d) phase 
recovery with modified HIO method with Huber penalty. 


application at hand. The results shown here have wide 
ranging applications in coherent imaging (optical. X-ray, 
electron microscopy, etc.) where the Fourier transform 
magnitude or the far-held diffraction intensity is read¬ 
ily available in experiments and the possibility of non- 
interferometric phase retrieval from intensity data can 
give access to much more valuable information about the 
object of interest. 
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